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Introduction 

Computational geometry is the branch of computer science 
concerned with the design and analysis of algorithms for geo- 
metric problems. Among the tools of computational geometry 
there seems to be a small set of techniques and structures that 
have such a wide range of applications that they deserve to 
be called fundamental, in the same sense that balanced binary 
trees and sorting are fundamental for combinatorial algorithms 
in general. In this paper we will discuss and illustrate a num- 
ber of these fundamental techniques of computational geometry. 
In most instances we give a worked out example of a problem 
solved by each of the techniques we discuss. 

Historical note. The term "computational geometry" was orig- 
inally coined by Robin Forrest [23] to denote the study of com- 
putational techniques in the realm of computer-aided geomet- 
ric design. More recently, this term has been used to name a 
somewhat different field, in a way broader and in a way nar- 
rower than Forrest's conception, a field which is now consid- 
ered part of theoretical computer science. This new use of the 
name "computational geometry," whose origin can be traced to 
Shamos [65], refers to the design and analysis of algorithms for 
all kinds of geometric problems, not necessarily in computer- 
aided design. This is the sense in which the term will be used 
in this paper. 

Influenced by other branches of theoretical computer sci- 
ence, this new field has focused on the design of efficient algo- 
rithms when the number of objects in the input is large. This 
consideration has led to the development of asymptotically effi- 
cient (and sometimes practical) algorithms for problems dealing 
with large numbers of simple underlying objects in two or three 
dimensions: points, lines, planes, polygons, and polyhedra. This 
situation is to be contrasted with that of computer-aided design, 
where the underlying objects are more complex, say bicubic sur- 
face patches, but where asymptotic considerations tend to be 
less dominant. Abstracting away from this dichotomy, we can 
say that the goal of computational geometry is to collect and 
study all techniques relevant to the computer description and 
manipulation of geometric objects, within the wider framework 
of the analysis of algorithms. 

It has been known since the time of Descartes that any geo- 
metrical problem can be recast in purely algebraic terms. It may 
therefore seem unnecessary to have a special discipline for geo- 
metric problems, distinct from numerical algebra and analysis. 
Indeed, some problems in geometry are best solved by algebraic 
methods, like computing the area of a polygon; but it is equally 
true that most geometrical concepts and algorithms are best 
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"seen" and studied in their own framework. Furthermore, most 
problems in computational geometry are neither purely combi- 
natorial, nor purely numerical, but rather an intimate mixture 
of both. Computational geometry has therefore its unique fla- 
vor, and its unique combination of tools and techniques. 

Many fundamental problems and results of computational 
geometry had been studied well before the field was recognized 
as a discipline in itself. Actually, we can say that computational 
geometry is the most ancient branch of computer science: the 
constructions of Euclidean geometry are legitimate algorithms, 
based on a well-defined, finite set of elementary operations. 

In fact, Euclidean geometry is where several of the key con- 
cepts of computer science were first introduced: the close rela- 
tionship between an algorithm and the proof of its correctness, 
the first examples of "provably unsolvable" problems (doubling 
the cube, trisecting the angle, squaring the circle), and so forth. 
Early this century the French mathematician Lemoine intro- 
duced (without much success) the idea of "computational com- 
plexity" of a geometrical construction, by counting the number 
of elementary steps it required. Other 19th century geometers, 
like Mohr and Mascheroni, showed that the ruler and compass 
of classical geometry could be replaced by other sets of tools 
(compass alone, ruler and scale, ruler and fixed-aperture com- 
pass, and so on). Their work is a close parallel to the theorems 
establishing equivalence of power for different flavors of Turing 
machines by simulation. 

Because of the long pre-history of the field and its large 
number of applications, the fundamental results and techniques 
of computational geometry are widely scattered in publications 
that range from highly abstract mathematics texts to appli- 
cations-oriented journals. It is only very recently that we are 
starting to see journals devoted explicitly (at least in part) to 
computational geometry, such as Discrete and Computational 
Geometry and Algorithmica. 

Applications. Computational geometry (in both senses) is 
currently undergoing rapid growth. The field now has attained 
some mathematical depth and maturity, yet it encompasses 
many fundamental questions that are far from being fully re- 
solved. Numerous areas of application contribute also to the vi- 
tality of the field; the list below shows a few of the most common 
ones, and some geometrical problems that are characteristic of 
each application. 

• Computer graphics 

geometric sorting in hidden surface elimination 
polygon intersection in clipping 
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geometric neighbor computations in hit detection 

• Computer-aided design 

union and intersection of geometric objects 

• VLSI design 

union, intersection, and near-contact of rectangles 
wire routing 

• Computer-aided manufacturing and robot control 

determining obstacle-free paths 
automatic milling 

• Pattern recognition and computer vision 

fitting of geometrical objects to noisy data 
clustering algorithms 

• Statistics, operations research, and numerical analysis 

classification by nearest-point determination 
finite-element decomposition 

In addition, there is an extensive list of problems and appli- 
cations that, although non-geometrical in origin and nature, 
can be better visualized and solved by recasting them in ge- 
ometrical terms. For example, the geometrical equivalent of a 
database range-query problem (like "find all employees with 
salary greater than $10,000 and age between 30 and 40") is the 
selection, from a collection of points in n-space, of those whose 
projections fall inside a given m-dimensional box, m < n. 

Eventually the techniques of computational geometry are 
bound to find more and more uses in computer-aided design, 
robotics, vision, and other applied areas. At the same time, 
these applied areas can serve as a rich source of problems for 
those in the field with more theoretical inclinations. It is the 
authors' hope that enough cross-fertilization will occur that in 
future years a distinction between the conceptions for the fields 
that Forrest and Shamos proposed will not need to be drawn. 

Overview. The paper is divided into three parts. In the first 
part we deal mostly with algorithm design methods. We have 
restricted our attention to techniques that are especially effec- 
tive in the geometric domain, that make use in some essential 
way of the geometry of the problem we wish to solve. In the 
second part we present a few data structuring techniques that 
are not themselves geometric but which have repeatedly found 
applications in geometric problems. Finally, in the third part 
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we discuss a couple of mathematically sophisticated techniques 
for the analysis of geometric algorithms. 

We have not tried to be exhaustive, and indeed we have 
completely ignored many important results and entire sub-areas 
of computational geometry. The interested reader will find more 
systematic coverage in the books by Preparata and Shamos [58] 
and Mehlhorn [49, vol. Ill], as well as in the surveys by Lee 
and Preparata [42], Edelsbrunner [20], and Dobkin and Sou- 
vaine [18]. Notable omissions are most range-query type prob- 
lems and the data structures motivated by them, such as seg- 
ment trees [4] and interval trees [46]; a good introduction to 
this area is Overmars' thesis [55]. Another important new tech- 
nique we have omitted is the use of randomization in geometric 
algorithms, for which the reader is referred to papers by Clark- 
son [13, 14]. We have also ignored the practical issues that arise 
in the implementation of computational geometry algorithms, 
and in particular the handling of "degenerate" cases where the 
input objects are not in "general position" (as theoretical dis- 
cussions usually assume them to be). Nevertheless we hope that 
our sampling is representative enough that the reader comes 
away with a sense of the technical excitement that permeates 
the field today. 



Geometric Algorithms 



5 



Part A: Algorithm Design Techniques 

In this first part we discuss techniques for the design of ge- 
ometric algorithms. Of course, computational geometry makes 
wide use of standard algorithm design methods, such as divide- 
and-conquer, dynamic programming, and so forth. However, we 
will concentrate here on methods and paradigms that seem to 
be more particular to this field, that in some essential way uti- 
lize the geometry of the problem at hand. Still, in many cases, 
the key insight for an efficient algorithm is an observation that 
allows us to map a geometric problem to a purely combinatorial 
one. We will see several instances of this in what follows. 



Al. The locus approach 



A common technique in solving geometric problems is the 
so-called locus approach. A classic example is the post office 
problem: given n sites in the plane (the post offices) and a query 
point x, report the site closest to x. If the sites remain fixed over 
several queries, then it pays to subdivide the plane into regions, 
each consisting of all points closest to a particular site. This 
partitioning of the plane is the well-known Voronoi diagram of 
the sites [58]. See figure 1. Once we have the Voronoi diagram, 
a nearest neighbor query can be answered simply by doing a 
point location in the diagram. This is the essence of the locus 
approach: subdivide space into regions such that all points in 
the same region yield the same answer to the type of query we 
are interested in. Point-location in this subdivision can then be 
used to answer any specific query. A method for computing the 
Voronoi diagram is given in section A4. 

One drawback of the locus method is that sometimes the 
size of the partitioning structure that we need to compute is 
too large. For example, in three dimensions, the Voronoi dia- 
gram technique is significantly less appealing for the post of- 
fice problem than it is in the plane, simply because a three- 
dimensional Voronoi diagram can have size 0(n 2 ). In such a 
situation one can trade space for time by building the patition- 
ing structure only for subsets of the given objects, and then 
querying each patitioning structure in sequence. For example, 
if we break the sites up into y/n groups of size y/n each and 
compute the Voronoi diagram for each group, then the total 
storage in three dimensions goes down to 0(ny/n), while query 
time increases by a factor of y/n. 

The locus approach is widely used in computational geom- 
etry algorithms, and plays a significant role in several of the 
examples to be given in the following sections. 




Figure 1. 



6 



Geometric Algorithms 



A2. Geometric transformations 





z 






fST~ 





Figure 2. 



Geometric transformations can be a very powerful tool in 
the development of geometric algorithms. In general, a geomet- 
ric transformation is an isomorphism between two geometric 
spaces (or a geometric space and itself). The simplest kind of 
transformation is a map that preserves the type and the essen- 
tial geometric properties of the mapped objects, such as the 
affine maps (which preserve parallelism), the projective maps 
(which preserve collinearity), and inversions (which preserve 
circles). Transformations of this sort can be used to simplify 
geometric algorithms, both in theory and in practice; for exam- 
ple, with a suitable affine map we can reduce many problems on 
an arbitrary triangle to problems on the equilateral one. Sim- 
ilarly, a projective map can be used to reduce many problems 
on a general conic section to problems on a circle. For more 
sophisticated examples, we can cite the linear programming al- 
gorithm of Karmarkar [37] and the minimum-area spanning el- 
lipse algorithm of Post [56]. Both of these arrive at the answer 
by transforming the problem through a sequence of geometric 
maps, which eventually reduce the problem to a trivial one. 

Even more useful are the geometric transformations that 
change the nature of the mapped objects. A transformation 
of this type is an isomorphism between two geometrical struc- 
tures, or two parts of the same structure, relating in a surprising 
way two classes of objects which appear to have very little in 
common. The isomorphism must by definition extend to the 
predicates and operations defined on those objects, and this is 
precisely where the the power of the method lies. By mechani- 
cally replacing the objects, predicates, and operations by their 
images, we can with negligible effort transport problems and 
solutions from one domain to the other. This not only reduces 
the number of theorems that need to be proved and of subrou- 
tines that have to be written, but also allows us to fully exploit 
in the study of one domain any intuition we may have gained 
in the other. 

Delaunay diagram via convex hull. An example of this class 
of transformation is the "lifting map" 

A : (x,y)>-> (x, y, x 2 + y 2 ), 

which lifts each point on the x,y- plane onto the paraboloid of 
revolution z - x 2 + y 2 . See figure 2. This map can be used to 
transform circular queries in the plane to half-space queries in 
three-dimensional space [33]. The reason is that any four points 
A, B, C, and D in the plane are co-circular if and only if A(A), 
\{B), X(C), and \{D) are coplanar. Therefore the intersection 
of the paraboloid of revolution z = x 2 + y 2 with any plane 
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meeting it projects into a circle in the x, j/-plane. Moreover, the 
part of the paraboloid below the plane projects to the interior 
of the corresponding circle in the x, y- plane, and the part above 
the plane to the exterior of the circle. 

Consider now n sites in the plane and their lifted images 
under the map A. These images are points on a convex sur- 
face, and therefore they are the vertices of a convex polyhedron, 
which is obviously their convex hull. Now consider a downward- 
looking face / of this polyhedron: by the definition of convex 
hull, its supporting plane x passes through the vertices of that 
face, and below all other vertices of the polyhedron. It follows 
that the projection of / onto the x, y-plane is a convex polygon 
(generally a triangle) whose circumcircle contains none of the 
other sites in its interior. We conclude that the projection of 
the downward-looking faces of the polyhedron are the faces of 
the Delaunay diagram of the given sites [7,27,33,40,41]. See 
figure 3. Similarly, an upward-looking face of the polyhedron 
corresponds to a convex polygon with vertices on the sites and 
whose circumcircle encloses all the sites. These are of course the 
faces of the dual of the furthest-point Voronoi diagram [58] for 
our collection of sites. 

It follows that algorithms for computing the convex hull 
of n points in three-dimensions yield, under this lifting map, 
algorithms for computing Delaunay triangulations (and there- 
fore Voronoi diagrams) in the plane. Curiously, it took a long 
time for the relationship between these two problems to be- 
come widely known, in spite of the fact that both have been 
solved independently by essentially the same algorithms. For 
example, the divide-and-conquer algorithms for Delaunay di- 
agrams by Guibas and Stolfi [33] and Lee and Schachter [43] 
are essentially the Preparata-Hong algorithm [57] for comput- 
ing the convex hull of n points in three dimensions, restricted 
to compute the lower half of the hull. 

The effect of the map A can also be obtained (at somewhat 
greater computational cost) by lifting the points of the xy plane 
onto a sphere with inverse stereographic projection [7]. The gen- 
eral idea behind the lifting map A is replacing non-linear tests 
on the input data (in our case, distance comparisons) by a non- 
linear map into a space of higher dimension, followed by linear 
tests on this transformed data. With a suitable lifting map, for 
example, we can reduce the problem of finding the minimum 
spanning ellipse of n points to a minimization problem on the 
faces of a 5-dimensional convex polyhedron [56]. Further ap- 
plications of such geometric transforms are discussed in Kevin 
Brown's thesis [7]. 




Figure 3. 
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A3. Duality 

Duality is a particularly important example of a geometric 
transform. In general, duality manifests itself as a non-trivial 
isomorphism between a mathematical structure and itself. As 
with any geometric transform, this isomorphism must extend 
to predicates and operations (and therefore to theorems and al- 
gorithms) denned on the objects. Familiar examples are the du- 
ality between variables and inequalities in linear programming, 
the laws of DeMorgan in boolean algebra, the current-voltage 
duality of circuit theory, the Fourier and Laplace transforms 
of analysis, and the face- vertex duality in the theory of planar 
graphs. 

As one might expect, this powerful concept has not es- 
caped the attention of geometers. A natural duality between 
points and lines in the projective plane has been known since 
the eighteenth century [15]. One way to express this duality is 
to map the point with Cartesian coordinates (a, 6) to the line 
with equation ax + by + 1 = 0, and vice-versa. We also map 
the point at infinity in the direction (a, b) to the line through 
the origin with equation ax + by = 0. Geometrically, if d is the 
distance between a point p and the origin O = (0,0), the dual 
of p is the line r perpendicular to Op and passing at distance 
1/d of O, on the side opposite to p. See figure 4. We will 
denote the line which is dual to point p by p*, and the point 
dual to line r by r*. (This is not the only possible duality be- 
tween points and lines, but we will stick with this version in the 
current section.) It is not hard to show that the mapping thus 
defined preserves incidence: if point p lies on the line r, then 
point r* lies on line p* . It also preserves any "betweenness" and 
continuity relations: as point u moves continuously from p to q 
along the line r, the line u* turns continuously at the point r* 
from p* to q*. 

Application to polygon visibility. For a typical use of du- 
ality in the development of geometric algorithms, consider the 
following visibility problem. Let P = (pi , P2 , ... p„) be a simple 
polygon in the plane, and let e denote the edge p\P2- Given a 
line r that enters P through the edge e, we wish to determine 
the first edge /(r) of P which r encounters after e. For con- 
creteness, we can imagine that P is a room with opaque walls, 
e is a luminous neon bulb, and we want to know which wall is 
illuminated by a given light ray emanating from e. See figure 5. 

Our approach to this problem is based on the observation 
that lines, and not points, are the primitive objects to con- 
sider in visibility questions. Since points are intuitively easier 
to grasp than lines, such questions are advantageously recast 
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in dual space, where the roles of points and lines are inter- 
changed. By using duality, we can obtain an algorithm which, 
after 0(n log n) preprocessing, allows us to solve this ray shoot- 
ing problem in 0(n) space and O(logn) response time. These 
latter bounds are clearly optimal. 

To simplify the discussion, let's assume P has been trans- 
lated, rotated, and scaled so that pi is the origin and pj is the 
unit point on the x-axis. Under the duality we defined above, 
pi is the line at infinity, and p\ 1S tne vertical line with equation 
x = -1. Now consider the set E of all lines intersecting e. It is 
easy to check that these lines dualize to all points with x < — 1, 
that is, the half-plane E* delimited by the line p| that does not 
contain the origin. See figure 6. 

Note again that the lines in E are in one-to-one correspon- 
dence with the light rays starting at e and directed into the 
polygon. Let us then classify each point p of E* according to 
the edge /(p*) that the ray p* illuminates. This defines a par- 
tition S(E*) of the half-plane E* into subregions where /(p*) 
is constant. Now consider two points p, q in the same region of 
S(E*), that is, with /(p*) = f(q*). This means the rays p* and 
q* illuminate the same edge of P. Let u be any point on the seg- 
ment pq; by the properties of duality, the line u* must intersect 
e, must be concurrent with p* and q*, and must lie in the angle 
delimited by p* and q* that does not include the origin. Since 
the boundary of P is connected and free from self-intersections, 
we must have /(«*) = /(p*) = /(?*); that is, u lies in the same 
region of S(E*) as p and q. We conclude the following: 

Lemma 1. Each region of the subdivision S(E*) is a convex 
polygon. 

Clearly the subdivision S(E*) (which has at most n — 1 
regions) provides the essential information we need to know in 
order to solve our problem. In order to compute this subdivision 
we first triangulate our polygon P by a standard 0(n log n) time 
algorithm [25]. We then consider each of the n — 2 triangles in 
turn, beginning with the one adjacent to the luminous edge e, 
and proceeding at each step across a diagonal from a visited 
triangle to a new one. In a typical step we enter a new triangle 
A = PiPjPk which has one side in common, say PiPj, with a 
previously handled triangle. At this point we know the region 
V of E* corresponding to all rays reaching p<pj (a diagonal of 
P) from e; this region is a convex polygon in the dual plane. 
These rays will now be subdivided into two groups, according 
to which other side of A they exit from. Dually, the line p£ cuts 
the polygon V into two subpolygons V',V", which are queued 
for processing in later steps. Assuming V is represented as a 
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balanced tree structure, we can perform each splitting operation 
in O(logn) time. So we have shown the following result: 

Lemma 2. It is possible to compute the convex subdivision 
S(E*) associated with illuminating the simple polygon P 
from edge e in O(nlogn) time and 0(n) space. 

Once we have the subdivision, standard point-location meth- 
ods can be used to solve the ray-shooting query in linear space 
and logarithmic time. 

Recently Tarjan and Van Wyk [72] have discovered that 
a triangulation can be computed in O(nloglogn) time. Also, 
by using a fancier representation for the polygons V that ex- 
ploits the finger search trees discussed in section B2, the split- 
ting operation described above can actually be carried out in 
linear total time. Thus the subdivision S(E*) can be built in 
O(nloglogn) time. 

A4. Space sweep techniques 

Space sweep is one of the most useful paradigms of computa- 
tional geometry. Generally speaking, space sweep is a technique 
that allows us to reduce an n-dimensional static problem to an 
(n— l)-dimensional dynamic problem. The basic idea can be de- 
scribed as follows. Suppose we have one or more objects lying 
in some Cartesian space, and imagine that space being com- 
pletely traversed by a moving hyperplane (the sweep plane). At 
any given moment, the sweep plane intersects some subset of 
the elements, the active ones. These intersections evolve contin- 
uously in time, except when certain discrete events occur: when 
new objects join the active set, when old objects leave it, or 
when the configuration of the intersections on the sweep plane 
undergoes a qualitative change. 

A space sweep algorithm is a discrete simulation of this 
process, using essentially two data structures: a queue of fu- 
ture events, and a representation of the active set and its cross- 
section by the current sweep plane. Each iteration of the algo- 
rithm removes the next^ event from the queue, and updates the 
cross-section data structure to mirror the effect of the sweep 
plane advancing past that point. Depending on the algorithm, 
each iteration may also reveal some future events that were not 
known at the beginning of the sweep; these must be inserted 
into the event queue, in the appropriate order, before proceed- 
ing to the next iteration. 
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The space sweep technique is quite old, and has probably 
been re-discovered several times in the early history of compu- 
tational geometry. It has been used in algorithms for the fol- 
lowing: finding planar convex hulls [58]; decomposing polygons 
into triangles [25], trapezoids [2], and monotone polygons [25]; 
intersecting line segments [5,44,45,66], polygons [5], and con- 
vex polyhedra [35]; merging convex maps [32, 44, 54]; computing 
Voronoi diagrams [24,27]; and many more. 

Voronoi via cones. We will describe here an optimal space 
sweep algorithm for computing the Voronoi diagram of n sites 
in the plane, due to S. Fortune [24]. His algorithm can be un- 
derstood more easily with the help of the following construction 
for the Voronoi diagram. For simplicity, we assume the sites are 
in general position, so that there are no two sites with the same 
x or y coordinate, and no four sites are cocircular. 

Imagine that an identical opaque cone with vertical axis 
is grown upwards from each site. See figure 7. A symmetry 
argument shows that any two cones intersect along a curve (a 
hyperbola) that is contained in a vertical plane. That plane 
is simply the perpendicular bisector in three-space of the two 
corresponding sites. 

Now imagine that we look at the collection of cones from 
below. The intersection of two cones will look like a straight 
line, the perpendicular bisecting line of the two sites p, q on the 
xy plane. Not every point of this line will be visible, however: 
a point a on it will be hidden by some other cone if and only 
if the site from which the latter grew is closer to a than p and 
q are. It follows that projecting the visible intersections down 
onto the plane gives the Voronoi diagram of the n sites. 

Sweeping the cones. Consider what happens when we sweep 
this three dimensional picture from left to right with a plane 
that is inclined at the same angle as the cones, with its top half 
trailing its bottom half. See figure 8. At any moment, the sweep 
plane meets the xy plane along a line which we call the sweep 
line. The intersection of a cone and the sweep plane is empty 
until the sweep line reaches the site p that is the apex of the 
cone. At that moment the sweep plane becomes tangent to the 
cone, and their intersection appears in the form of a degenerate 
parabola of zero width, whose projection is a ray starting at p 
and pointing towards x = — oo. A moment later that ray opens 
up to become a narrow parabola, which grows wider and wider 
as the sweep plane moves further to the right. It is easy to see 
that the projected parabola on the xy plane will have the site 
p at the focus and the sweep line as the directrix. Referring to 
figure 8, we see that triangles pvw and uvw have equal angles 
and share an edge; therefore they are congruent and pw = uw. 




Figure 7. 
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The parabolic front. Now let's make the sweep plane opaque, 
too, and look at the whole thing from below. The sweep plane 
will completely hide the cones whose sites lie to the right of the 
sweep line. The cones of the sites already swept over will give 
rise to a set of left-pointing parabolas with parallel axes and 
various widths. A point on the xy plane is inside (to the left 
of) a parabola if and only if the cone of the latter lies below 
the sweep plane at those (x,y) coordinates. Therefore, any part 
of a parabola that lies inside another one is invisible. The only 
visible parts of those parabolas are those on the boundary of 
the union of their interiors. We call this boundary the parabolic 
front. See figure 9. 

Note that any horizontal line intersects the parabolic front 
at exactly one point. Therefore, the front consists of a single 
chain of arcs, extending from y = — oo toy = +oo. To the right 
of the parabolic front, the sweep plane lies below all the cones, 
and completely blocks the view from the underside. To the left 
of the front, the plane lies above one or more cones, and cannot 
hide anything that would be visible otherwise. Looking from 
below, what we see there is a cut-out portion of the Voronoi 
diagram of the sites. It is easy to see that each arc of the front 
lies in some Voronoi region, and each "break" between two con- 
secutive arcs lies on a Voronoi edge. 

Evolution of the parabolic front. As the sweep plane ad- 
vances, the arcs in the parabolic front move to the right and 
become flatter, and the breakpoints between them move along 
the corresponding Voronoi edges. It is not hard to show that 
the more recent of two consecutive arcs (the one whose site lies 
further to the right) grows in vertical extent at the expense of 
its older neighbor. We claim that as the sweep plane moves from 
x = — oo to x = -foo, the breakpoints of the front will trace out 
every edge of the Voronoi diagram: 

Lemma 3. Every point of every edge of the Voronoi diagram is 
a breakpoint of the parabolic front at some time during the 
scan. 

Proof: Let e be any Voronoi edge, and u a point on e (other 
than its endpoints). Consider the situation when the sweep line 
is already to the right of u, and the distance r between the two is 
equal to the distance between u and the sites p, q of the Voronoi 
regions separated by e. See figure 10. From the properties of 
Voronoi diagrams we know that no other site is closer to u than 
those two, that is, there are no other sites on or within the circle 
C with center u and radius r. Obviously, at that moment u will 
be on the intersection between the parabolas associated with p 
and q. 
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The only way u could not be a breakpoint of the front is 
if there were some other parabola passing between u and the 
sweep line. Let v be the point on that parabola just to the right 
of u; the focus of the parabola must then be on the circle with 
center v that is tangent to the sweep line. But that circle lies 
inside the circle C, a contradiction. We conclude that u is on 
the current front. | 

The discrete events. The continuous evolution of the front is 
punctuated by discrete events, where new arcs join the front, or 
old arcs shrink down to nothing and drop out of the picture. The 
first case happens, for instance, when the sweep line runs over 
a new site s. See figure 11. That introduces a new degenerate 
parabola a with zero width (i.e. a twice-traversed horizontal 
ray) with focus and apex at s. The portion of o that lies to the 
right of the parabolic front becomes a new arc of the latter, and 
its insertion will split an existing arc (3 of the front in two pieces 
(if the sites are in general position, we can ignore the possibility 
of the ray passing through a breakpoint). We call this a site 
event, and we make the following claim: 

Lemma 4. The only way for an arc to enter the parabolic front 
is through a site event. 

Proof: The only other alternative is for new arcs to arise due to 
changes in the shape and position of existing parabolas, that is, 
due to some parabola overtaking the front and breaking through 
it. However, this cannot happen. Consider the situation when 
the parabola a in question is about to break through. It cannot 
do so in the middle of an arc /3, since that would require the 
curvature of a to be higher than that of f3 and the focus of a 
to be further from the sweep line than that of f} — which are 
contradictory statements. 

On the other hand, a cannot break through the corner be- 
tween two arcs (3 and 7. To do so, it would have to move faster 
in the i-direction than (3 and 7 at the y-coordinate of their in- 
tersection. A little algebra shows that the horizontal speed at 
which a particular parabola advances at a given y-coordinate 
is a monotonically decreasing function of its absolute slope \a\ 
(more precisely, it is |(1 + 1/er 2 )). Therefore, for a to overtake 
/3 and 7 at a point u, the absolute slope of a at u would have 
to be smaller than those of (3 and 7, contradicting the fact that 
the latter are consecutive arcs of the front meeting at u. | 

Let's now turn to the second class of events, where an ex- 
isting arc drops out of the front after having shrunk down to a 
single point 11. At that moment there are three parabolas pass- 
ing through u, namely the disappearing arc f3 and the arcs a, 7 
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just above and below it on the front. See figure 12. As in the 
proof of lemma 4, it is easy to show that a and 7 cannot belong 
to the same parabola. The point u is therefore equidistant from 
the three corresponding sites and from the current sweep line. 
In other words, the sweep line is tangent to (and to the right 
of) the circumcircle of the three sites. We conclude that the 
point u is a vertex of the Voronoi diagram. (Since the sites are 
in general position, we can ignore the possibility of four or more 
of them having the same circumcircle.) We call the rightmost 
point w of the circumcircle associated with a Voronoi vertex u 
a circle event. From the preceding discussion, the following 
result is clear: 

Lemma 5. The only way for an arc to leave the parabolic front 
is through a circle event. 

The data structures. Now that we understand how the para- 
bolic front evolves during the sweep, we are ready to consider its 
simulation in the computer. As mentioned before, the algorithm 
needs two main data structures: a queue of future events, and 
a description of the current parabolic front. The former con- 
tains all site events that lie to the right of the sweep line. It 
also contains some of the circle events in that region, namely 
those that correspond to three consecutive arcs of the current 
parabolic front. Note that some future circle events will be 
missing from this list, either because they involve sites that 
haven't been swept over yet, or because the corresponding arcs 
are not yet consecutive elements of the parabolic front. As we 
will see, these circle events will be discovered during the sim- 
ulation, and inserted in the event queue before the sweep line 
advances past them. Conversely, not every three consecutive 
arcs of the current front specify a circle event. 

The parabolic front is represented by a balanced search tree 
(or any equivalent structure), with the arcs as leaves and the 
breakpoints as internal nodes. We will use this tree to efficiently 
insert and delete arcs of the front, to locate the arc that inter- 
sects a given horizontal line, and to locate the arcs immediately 
above and below a given one. 

We also hang from this tree the part of the Voronoi data 
structure built so far, consisting of the vertices, faces, and edges 
that lie totally or partially to the left of the current parabolic 
front. Each leaf of the tree (i.e., each parabolic arc) includes a 
pointer to the corresponding site, and each internal node (break- 
point) includes a pointer to the record representing the associ- 
ated Voronoi edge in the incomplete diagram. 
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Event processing. Given these data structures, the algorithm 
itself is a straightforward event-driven simulation loop. At each 
iteration we remove the next event from the queue (the one with 
smallest x-coordinate), and we simulate the effect of the sweep 
line advancing to that point. We again use the assumption of 
the sites being in general position to ignore the possibility of 
two events having the same x-coordinate. 

To process a site event s, we add to the tree a new degen- 
erate arc a, consisting of a left-pointing horizontal ray starting 
at 3. This requires us to split an existing arc /? in two, and 
add two new breakpoints. We must also add a new edge to the 
Voronoi diagram, associated with both breakpoints, which is 
the frontier between the Voronoi regions of s and of the site 
associated with f3. See figure 11. 

Processing of a circle event is equally simple. Recall that 
such an event corresponds to the sweep line reaching the right- 
most point of the circumcircle of three sites a, b, c, associated 
with three consecutive arcs a, /?, 7 of the parabolic front. At 
that point the middle arc f3 has zero length and must be deleted 
from the tree. Before we do so, we must add to the Voronoi dia- 
gram a new vertex u, the circumcenter of the triangle abc. That 
vertex is incident to the two Voronoi edges currently associated 
with the breakpoints a-/3 and (3-^. Also, after deleting /? we 
must add a new Voronoi edge incident to u, and attach it to 
the new breakpoint a-7. See figure 12. 

Event scheduling. We still haven't said how and when the 
events get inserted into the queue. Site events get inserted all 
at once, at the beginning of the algorithm. Circle events are 
inserted as the simulation proceeds, when the three associated 
arcs first become consecutive elements of the front. 

More precisely, when adding a new arc /? during a site 
event s we look at the two new consecutive arc triplets created 
by the insertion (having (3 as the first and last element, re- 
spectively). For each triplet we compute the circumcircle of the 
corresponding sites, and add its rightmost point to the queue 
as a circle event. Note that since the site s is on the current 
sweep line, the new circle event cannot lie to the left of that 
line. 

Similarly, when deleting an arc /? between arcs a and 7 
during the processing of a circle event, we check the two new 
consecutive arc triplets created by the deletion (which start at 
a and end at 7, respectively). For each triplet we compute the 
circumcircle of the sites, and add its rightmost point to the 
queue, but only if it lies to the right of the sweep line. Note 
that it is possible for the circumcircle to lie strictly to the left 
of the sweep line, as figure 13 shows. That means the center 
of the circumcircle has already been swept over by the front. 




Figure 13. 
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As discussed below, that center (if it was a Voronoi vertex) will 
have been noticed and processed at an earlier date. 

It is easy to see why this algorithm schedules every discrete 
event affecting the front in due time, that is, before the sweep 
line reaches that event. Because of lemma 4, new arcs can enter 
the front only through site events, which are all scheduled 
beforehand. Because of lemma 5, whenever an arc is about to 
disappear from the front, the sweep line must be tangent to the 
circumcircle of the sites associated with it and the two adjacent 
arcs. Therefore, the three arcs must have become adjacent at 
some earlier time, and the corresponding circle event must 
have been inserted in the queue at that moment. 

On the other hand, the algorithm may schedule spurious 
circle events which do not correspond to actual changes in 
the front, and therefore to actual Voronoi vertices. This may 
happen, for example, if three arcs become consecutive elements 
of the front at some point, but a new arc is inserted in their 
midst before the sweep line reaches their circle event. At the 
moment the circle event was inserted, the Voronoi edges be- 
tween the three arcs looked as if they would converge at some 
later point, but the arrival of the new site caused those edges 
to terminate prematurely. See figure 14. In fact, by the time 
the circle event is reached, some of the three arcs may have 
already disappeared, engulfed by more distant neighbors or by 
new arcs inserted in their midst. 

To fix this problem, we keep for each parabolic arc /3 on the 
front a pointer death(/3) to an event in the queue, namely the 
earliest circle event where that arc is currently scheduled to 
disappear. The Voronoi vertex corresponding to that event is 
the predicted meeting point of the Voronoi edges being traced 
by the endpoints of /?. If those edges diverge, the arc is currently 
expected to live forever, and its event pointer is nil. To avoid 
spurious circle events, it suffices to keep the death pointers 
up to date during the simulation, and to promptly remove from 
the queue any circle event that is no longer pointed to by 
any arc. More precisely, we must discard the event death(a) 
whenever the arc a is split in two by a site event, or whenever 
one of the two arcs adjacent to a is deleted by a circle event. 
This completes the description of Fortune's algorithm. 

Analysis. In order to get the desired bounds on the space and 
time cost of the algorithm, we have to show that at any time 
the parabolic front and the event queue contain at most 0(n) 
elements. This is not entirely obvious, since, for example, each 
parabola may give rise to many arcs on the front. However, the 
only way the parabolic front increases in size is through site 
events: then a new arc is added and an old arc is split in two. 
Since there are only n such events, this proves that the front 
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has linear size. (Another way to derive linearity comes from the 
fact that the front cannot contain two parabolas a and {3 that 
each appear twice in the pattern ... a .../?... a .../?.. .. The 
Davenport-Schinzel sequence theory developed in section C2 
then implies that the combinatorial complexity of the front is 
0(n).) Since each event currently in the queue is either a site 
event or a death event for an arc on the front, we conclude that 
the queue too never contains more than 0(n) elements. With 
a suitable priority queue structure [1] we can insert and delete 
an event in the queue, or find the one with minimum x, in only 
O(logra) time. We conclude that each iteration takes O(logn) 
time, and therefore we arrive at the following result: 

Theorem 6. Fortune's algorithm computes the Voronoi diagram 
of n sites in 0(n log n) time. 



A5. The configuration space approach 

The configuration space approach is a geometric technique 
from robotics research that has found application in many other 
areas of computational geometry [63]. In a typical application, 
we want to check quickly whether two or more geometric objects 
(say, a robot and the furniture in a room) intersect when placed 
in a given relative position. If we have to perform this test for 
the same set of objects in many different positions, we may con- 
sider pre-computing the entire set 5 of all relative placements 
(configurations) for which the two objects intersect. Then each 
query reduces to a test whether the proposed placement lies 
within S, that is, to a point location problem in the space of all 
configurations. 

This approach is practical only when the configuration space 
has low dimension, such as problems involving two rigid objects 
in the plane. As the objects' motion becomes less constrained 
(allowing, say, translations and rotations in three-space, or ar- 
ticulated joints) the dimensionality of the configuration space 
increases, and the cost of pre-computing, storing, and searching 
S grows exponentially. 

The simplest case that leads to non-trivial algorithms is 
detecting intersections between two polygonal objects P and Q 
in the plane that are allowed to move only by parallel trans- 
lation, without rotation. In this case the configuration space is 
two-dimensional, and consists of all displacement vectors mea- 
sured between two reference points fixed on the objects. The set 
S of configurations leading to "interference" between the two 
polygons is itself a planar polygonal region. 
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Even with these restrictions, the proper treatment of this 
problem for general polygons requires more machinery than we 
can afford to develop here. The reader is referred to the litera- 
ture, in particular to the kinetic framework proposed by Guibas, 
Ramshaw, and Stolfi. [31]. Henceforth we will restrict ourselves 
to the case of two convex polygons. 

The set of all invalid placements can be easily expressed 
in terms of the Minkowski sum or convolution of two sets A, B 
of vectors, which is simply the set of all pairwise sums of their 
elements: 

A* B = {a + b : a G A and b G B} 

To see how this relates to our problem, we must consider the 
two given objects as sets of vectors (measured from a common 
origin). Let's also denote by A N the set A rotated 180° around 
the origin, and by A x the set A translated by the vector x: 

A N = {-a : a G A} 
A* = {a + x : a G A} 

The set of forbidden positions is given by the following trivial 
result: 

Lemma 7. If P and Q are plane objects, then P x intersects 
if and only if the vector x — y lies in the set P N * Q (or, 
equivalently, y — x lies in P * Q^). 

In other words, P N * Q is the set of displacements of P relative 
to Q for which the two figures intersect. 

The convolution of two polygonal figures P, Q is also a 
polygonal figure. If P and Q have p and q edges, respectively, 
then P * Q (and P N * Q) may have Q(pq) edges in the worst 
case. As figure 15 shows, P * Q need not be simply connected, 
even when P and Q are simple polygons. If P and Q are convex, 
however, P * Q is a convex polygon with at most p + q edges. 

It is convenient to view each edge of a polygon as a vector 
directed counterclockwise around the polygon. Observe that the 
edges around a convex polygon are ordered by "slope" (direc- 
tion angle). In the case of convex polygons, P * Q has a simple 
characterization: its edges are those of P and Q, merged in slope 
order. More precisely, if we imagine the edges directed counter- 
clockwise, then every edge e of P is translated by the vertex v 
of Q such that the direction of e lies in the exterior angle at v 
that is facing counterclockwise. See figure 16. 

This fact allows us to compute the convolution P N * Q 
using only 0(p + q) time and space. After that we can test in 
0(log(p + q)) time whether P displaced by a given vector x 
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would intersect Q displaced by y, by simply checking whether 
x - y is inside the convex polygon P N * Q. We conclude the 
following: 



Theorem 8. It is possible to pre-process two convex polygons P 
and Q of size n using only O(n) time and space, so that any 
translated copies of P and Q can be tested for intersection 
in O(logn) additional time. 



Intersection without preprocessing. We may be tempted 
to rest on our laurels at this point. After all, don't we have to 
look at each edge of P and Q to know if they intersect? Surpris- 
ingly, the answer is no. By a suitable kind of binary search, we 
can test for intersection in time only 0(\og(p + q)), without any 
pre-processing. Instead of precomputing the convolution P™ * Q 
explicitly, we compute each edge of the convolution only if and 
when it is needed. 

Imagine that we move P horizontally to the left, all the way 
to infinity. We call the region swept over this way (including P 
itself) the left shadow of P, denoted by Pl- See figure 17. The 
right shadow Pr or P is similarly defined by sweeping P to the 
right. Obviously, P = P L n Pr. It is easy to verify the following 
result: 

Lemma 9. The intersection P C\Q is non-empty if and only if 
both Pl n Qr and Pr fl Ql are non-empty. 

Thus, to check whether P x intersects Q v we need only know 
how to check for intersections between a left shadow L u and a 
right shadow R v , that is, between {Pl)* and (Qr)*, or between 
or {Ql) v and (Pr) x . By lemma 7, this is equivalent to testing 
whether v — u lies in L * i2 N . Note that i2 N is a left shadow, and 
that the convolution of two left shadows is also a left shadow. 
This is then the problem we will consider: given two left shadows 
A and B, discriminate a given point w against A * B. 

As stated earlier, the edges of the convolution A*B (viewed 
as vectors directed counterclockwise around the figure) are pre- 
cisely those of A and B, merged in order of direction. Suppose 
the direction of an edge e from A lies between the directions 
of two consecutive edges f,g of B, both incident to the same 
vertex v. Then the edge e appears in the convolution A * B 
displaced by the vector v. In the same way, each edge of B gets 
displaced by some vertex of A. 

In what follows, placing e in the convolution means deter- 
mining this vertex v, and hence the position of e in A * B. 
Let ua and be the number of edges in A and B, respec- 
tively, and let n = + be the number of edges in the 
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convolution. If the coordinates of the vertices of each chain are 
stored in a linear array, in counterclockwise order, we can use 
binary search to place any edge e in the convolution, in time 
(^(maxllogn^lognB}) = O(logn). Once e has been placed, 
we can test w against its supporting line / and its endpoints 
p,q in 0(1) additional time. See figure 18. If w lies to the 
right of /, then w is outside A* B, and we can stop. If w lies in 
the left shadow of e, then it definitely lies inside A* B, and we 
can stop. 

If neither case holds, then we must continue checking w 
against other edges of A * B; specifically, either the edges that 
follow e or those that precede e in counter-clockwise order, de- 
pending on whether w lies above or below e's shadow. Suppose 
e is an edge from A that was displaced by a vertex v of B, and 
w lies above its shadow; then the edges of A * B we must con- 
sider next are those that follow e in A, plus those that follow 
v in B. The other cases are symmetrical. In other words, we 
have managed to transform the original problem into a simi- 
lar (but smaller) one: discriminate w against the convolution of 
two convex left shadows, each defined by a string of consecutive 
edges from A or B. By repeating this process we will eventually 
reduce the convolution to a single edge, which either shadows 
w or leaves w to its right. 

How many iterations will this take? If the edge e we dis- 
criminate against is always the median edge of A * B, then at 
each iteration the problem size gets halved, and the number of 
iterations will be at most log 2 n. Finding the median edge of 
A * B is a bit tricky, given that A and B are stored in separate 
arrays. Nevertheless, by a binary search procedure it is possi- 
ble to find this median in 0(log n) time. We will not bother to 
give the details, however, since a much simpler method gives 
the same asymptotic result. 

Note that it is not necessary to find the exact median of 
A * B to guarantee a logarithmic number of steps. It suffices to 
choose e so that at each stage we always eliminate at least some 
fixed fraction a of the edges. In particular, if we let e be the 
median of the longest of the two chains (which can obviously be 
found i onstant time), then at least 1/4 of the edges of A * B 
must come before e, and at least 1/4 must come after it. There- 
fore at each stage the size of the problem is reduced by a factor 
< 3/4. The number of stages will be at most log 4/3 n, which is 
bigger than log 2 n by a factor of log(4/3)/log2 = 2.409+, but 
is still O(logn). Since each stage (placing e in the convolution 
and testing w against it) takes O(logn) time, it follows that 
testing whether w lies in A * B can be done in O((log ra) 2 ) time. 
We conclude the following: 
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Theorem 10. With no preprocessing, it is possible to test whether 
two convex polygons of size n intersect in O((logn) 2 ) time. 

Intersection in O(log n) time. Surprisingly, we can still im- 
prove on this result. We will show that we can actually discrim- 
inate the point w against the convolution A * B in time only 
O(logn). To do this, we cannot afford to spend time O(logn) 
on every test just to place the edge e in the convolution. Instead 
we wish to find a constant time test that allows us to get rid of 
some fraction of the edges of A or B. With this in mind, let us 
look at the convolution from yet another viewpoint. 

Let a', a" be the lower and upper endpoints of the chain 
A, and let b',b" be those of B. The convolution A* B will be 
a right-convex path on the plane from a' + b' to a" + b". Now 
consider all possible ways to interleave the edges of A with the 
edges of B, while maintaining the relative order of any two 
edges from the same chain. Each interleaving defines a path on 
the plane from a' + b' to a" + b" (since all paths consist of the 
same set of vectors, their total displacement is the same). Every 
path consistently moves in the general direction of increasing y, 
but may turn left and right several times. Only one path will 
be convex to the right: the one in which all edges are sorted by 
direction, that is, A* B. 

Note that A * B is also the rightmost of those paths, in the 
sense that no point on any other path lies to the right of A * B, 
and any other path has at least one point that is strictly to the 
left of A * B. To see why this is true, observe that if / and g are 
consecutive edges in a path H and are not in the right order, 
then swapping them will produce another valid path H' that 
lies to the right of H. See figure 19. 

Now let / denote the median edge of A, and g the median 
edge of B. Without loss of generality, we can assume that / 
precedes g in slope order, and therefore also in the convolution 
A * B. Let Ah, Al, Bfj, Bl denote respectively the high and 
low halves of A and B that are separated by / and g. 

Consider next the chain D which consists of Al * Bl, fol- 
lowed by /, and g, followed Ah*Bjj, in this order. See figure 20. 
Since the displacement of Al * Bl is the sum of the displace- 
ments of Al and of Bl, the position of / in the chain D can be 
computed in constant time. 

Let us then test w against the edges / and g, positioned 
as in the path D. If the point w falls in region LEFT (see 
figure), then clearly w is to the left of the convolution, and we 
are finished. If w falls in region ABOVE, then what can we 
conclude? Observe that in the convolution A * B, the edge / 
and all the edges in Al must lie below the dashed line. The 
reason is that / precedes g in slope order, hence it precedes 
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also all edges of An and Bfj. Therefore, the path leading to 
/ in the convolution is a subset of Al * Bl, which is the path 
leading to / in D. 

Therefore, if w falls in region ABOVE, we can discard / and 
all edges in Al, since they cannot affect the classification of w. 
Similarly, if w falls in region BELOW, then in the convolution 
w will be below g and all edges of Bh, so these edges can be 
discarded. By this constant time test we are able to eliminate at 
least half of the chain A or half of the chain B. We can repeat 
this step until one of the chains (say, A) gets reduced to a single 
vertex v. 

Note that we can eliminate only half of one chain, and 
we don't get to choose which one. Since the two chains may 
have very different lengths, we cannot guarantee that the size 
+ n B of the convolution will decrease by a constant fraction 
at each stage. However, after 0(log tia) + O (log ns) = O(log n) 
such steps one of the two chains (say, A) will be reduced to a 
single vertex v. This means the convex shadow bounded by that 
chain is a single horizontal ray extending from v to the left. The 
convolution of the two chains is therefore the other chain (B) 
displaced by the vector v. Checking w against the convolution 
is then equivalent to checking w-v against the chain B, which 
can be done in O(logn) additional time by a straightforward 
binary search. The total time for checking w against A * B is 
therefore 0(log n) +0(log n) = O(log n). So we have proved the 
following result: 

Theorem 11. It is possible to test if two convex polygons of size 
0(n) intersect in time O(logn). 

The same result was obtained by Chazelle and Dobkin, us- 
ing different methods [9]. 

To conclude this example, let's comment briefly on the form 
in which the answer should be presented. In many practical 
applications, just knowing whether P and Q intersect is not 
enough. We often need a witness or certificate that corroborates 
the answer: a common point if the polygons do intersect, or a 
separating line if they don't. Our algorithm for testing whether 
a point w lies in the intersection of two convex shadows A and 
B returns such a witness, namely an edge e of A*B such that w 
is either in the shadow of e or to the right of its supporting line. 
By that time we also know the chain where e came from, and 
the vertex v of the other chain by which e was displaced. From 
these data we can easily recover a witness for the intersection 
or disjointness of the original polygons. 
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A6. Separator techniques 



As in graph theory, the existence of "small" separators often 
leads to efficient divide-and-conquer algorithms. Perhaps the 
best known result of this kind in computational geometry is 
Chazelle's polygon cutting theorem [8]. That theorem states 
that a simple polygon admits of a diagonal that cuts it into two 
subpolygons such that (roughly) the larger is not more than 
twice as big as the smaller. See figure 21. (As usual, the size 
of a polygon is its vertex count, not its length or area.) This 
kind of balanced decomposition is crucial when we seek efficient 
divide-and-conquer methods. 

The existence of such a diagonal for a simple polygon fol- 
lows from a corresponding theorem for trees. Triangulate the 
polygon and look at the tree T that is the dual graph of the tri- 
angulation. The nodes of T have degree three; by a well-known 
result of graph theory, T has an edge whose removal cuts it in a 
balanced fashion. The dual of this edge is the diagonal we want. 
Chazelle showed also that such a cutting diagonal can be found 
in time proportional to the size of the polygon. 

Suppose that our polygon P has n sides. If we cut P along 
a cutting diagonal and then recursively apply the technique to 
the two subpolygons thus formed, we will obtain a tree S which 
forms a balanced hierarchical decomposition of P. The leaves of 
this tree are the triangles of a triangulation of P, while internal 
nodes correspond to diagonals used in the triangulation. The 
whole process can be carried out in O(nlogn) time. If all the 
produced subpolygons are saved within the appropriate node of 
5, then the total storage will also be 0(n log n). Such a balanced 
hierarchical decomposition is very useful in computing shortest 
paths inside P [29], or in visibility questions inside P [11]. In 
most applications the redundancy involved in storing all the 
subpolygons can be removed with some cleverness so the storage 
used can be made linear in n. 




Figure 21. 



A7. Decimation 

Decimation is our name for a programming technique due 
to Nimrod Megiddo that has yielded surprisingly efficient al- 
gorithms for a variety of computational geometry and opera- 
tions research problems. The technique has also been called the 
Megiddo method or prune-and-search by other authors. Some 
applications are computing the smallest enclosing circle [47] and 
the convex hull [39] of n points on the plane, and solving linear 
programming problems in spaces of fixed dimension [48]. 

Problems where the decimation technique is applicable seem 
to be those where the answer is some simple geometric object 
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Figure 22. 






Figure 23. 



that ultimately depends on only a few of the input data. The 
technique depends on the existence of some quick criterion that 
allows us to eliminate in linear time some fixed fraction of the 
input data. By repeatedly applying this process with the re- 
maining data, we will be left eventually with a trivial problem 
that can be solved efficiently by some other method. If each 
pass reduces the size of the problem by a factor p < 1, the total 
time required by the algorithm will be 0(n + pn + p 2 n + •■•), 
which is still O(n). 

A good example is that of computing the smallest circle 
enclosing n given points in the plane. Such a circle is always 
determined by either two or three of the given points. This 
classic problem of computational geometry was conjectured to 
have ft(nlogn) worst-case cost, until Megiddo found a decima- 
tion method that on each pass eliminates a fixed fraction of the 
points which do not support the minimum circle [47]. 

The common tangent problem. We will describe here in 
detail another and somewhat easier example, the computation 
of the common upper tangent of two sets of points. The prob- 
lem is this: we are given a collection P of n points in the plane, 
which is separated into two non-empty subsets L and R by a 
known vertical line m, with L on the left and R on the right. We 
wish to find a line t passing through one point from each subset, 
such that none of the given points lies above t. See figure 22. 
In other words, we want the upper exterior common tangent of 
the convex hulls of L and R. 

If the convex hulls of L and R were known, the common 
tangents could easily be found in linear time [57]. However, com- 
puting the convex hull of n points costs 0(nlogn) in the worst 
case [73]. Is this also a lower bound to computing the common 
tangent? The answer is no: as Kirkpatrick and Seidel showed, 
the upper common tangent of L and R can be computed in 
time linear in n, without knowledge of their convex hulls [39]. 
The key to their algorithm is a simple decimation criterion that 
in linear time allows us to eliminate a good many of the points 
that do not define the common tangent. 

Let us fix for the moment our attention on lines of a par- 
ticular slope a. We can compute in linear time a point p of L 
that is extremal for this slope, i.e. such that a line through p 
with slope a leaves all other points of L below it. We can do the 
same for R and obtain a corresponding extremal point q. Now 
if the line pq has slope less than a, then so must the common 
tangent t; similarly, if pq has slope greater than a then so does 
t; and if pq has slope equal to a then t = pq. See figure 23. 

Now suppose the first case holds, so the slope of t is less 
than a. Let r,s be any two points of L U R, such that r is to 
the left of s and the line rs has slope greater than a. Then we 
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can conclude that t cannot pass through r, because a line of 
slope less than a through r must pass below s. See figure 24. 
The second case, where the slope of t is greater a, is entirely 
symmetrical: we can eliminate the second member of any pair 
r,s, with r to the left of s, if the line rs has slope less than a. 

These remarks suggest the following decimation method. 
Pair up the n given points in an arbitrary way, and find the 
median slope a of the n/2 lines defined by those pairs. Now 
compute the extremal points p and q for the slope a, and elim- 
inate one point from every pair that satisfies the criteria de- 
scribed above. It should be obvious why the median is a good 
choice: since half the pairs have slope less than a, and half have 
slope greater than a, then half the pairs will satisfy the crite- 
rion, no matter whether the slope of pq is greater or less than 
a. So, in either case we eliminate one point from half the pairs. 
Of course if pq has slope a then we are done. 

It is possible to find the median in time linear in n using the 
rather elaborate technique of Blum and others [6]. The other 
operations clearly take 0(n) time. Therefore, in linear time we 
either stop, or eliminate 1/4 of the original points. If we repeat- 
edly apply this elimination process on the remaining points, we 
are guaranteed to find the common tangent, at a total cost of 
0(n + (3/4)n + (3/4) 2 n + •••) = O(n) Note that we may elimi- 
nate a different number of points from L than from R, but this 
does not affect the analysis. 

Theorem 12. The upper exterior common tangent of two ver- 
tically separated point sets can be computed in linear time. 
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Kirkpatrick and Seidel used this result to obtain a planar 
convex hull algorithm that runs in time 0(n log h), where h is 
the size of the hull [39]. The one step that makes the algorithm 
above impractical is finding the median. Unfortunately, such 
a step is almost always present in decimation algorithms. For 
practical purposes we can substitute a randomized median (or 
other quartile) finding method. This yields algorithms that are 
linear-time only in a probabilistic sense but are much simpler 
to implement. 



Figure 24. 
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A8. Hierarchical structures 

Several geometric problems have been successfully solved 
by the following approach. Given an instance P of the problem, 
we may attempt to extract from it a "coarsened sample" P', 
another instance of the same problem such that (1) the size of 
P' is only a fraction of the size of P, and (2) a solution to P' 
yields a solution P with only little additional effort. If we can 
define and effectively construct such samples, we can solve P 
by sampling it, in turn sampling the sample, and so on, till we 
obtain an instance small enough to be solved by some trivial 
method. We then back up through the sample hierarchy we 
have constructed, obtaining the solution for increasingly finer 
samples, till we come to the original problem P. 

David Kirkpatrick was the first to popularize such hierar- 
chical decomposition methods in computational geometry. He 
used this approach in the first practical logarithmic point loca- 
tion method for planar subdivisions [38]. He and Dobkin also 
gave many applications to intersection questions about convex 
polyhedra [16, 17]. 

Figure 25. A segment visibility problem. Let's illustrate this hierar- 

chical decomposition technique on a very simple problem. Sup- 
pose we are given a set S of n vertical line segments in the 
plane. For ease of exposition we assume that the x-coordinates 
of all these segments are distinct, so no overlaps among them 
are possible. We are interested in organizing these segments into 
a structure that allows us to answer efficiently queries of this 
form: given a point p in the plane, report the segments in S 
(if any) immediately to its left and to its right. In other words, 
if we extend horizontal rays from p left and right, we wish to 
know the first segments of S that will be hit. See figure 25. 

Notice that this is a point location problem. From each 
segment endpoint, let's extend rays to the left and right till we 
encounter another segment of S. The segments of S together 
with these rays form a rectangular partition of the plane. See 
figure 26. All points lying in the same rectangle of this par- 
tition yield the same answer to our "left- and- right neighbor" 
query. By sweeping the plane with a horizontal line, we can con- 
Figure 26. struct this rectangular partition in O(relogn) time. We could 

then solve our problem by triangulating the regions and using 
Kirkpatrick's original point location method. It is preferable, 
however, to apply the sampling idea directly to the original 
problem, so as to avoid the complexity of the triangulation al- 
gorithm and the numerical difficulties involved in testing a point 
against oblique lines. 
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Sampling the segments. How are we to define a sample of 
57 Let us name the segments of 5 as s\ , s 2 , . . . , s n in some 
order. We say that segment s, "sees" segment Sj if there is some 
ordinate at which these two vertical segments are adjacent in 
the horizontal ordering. In more informal terms, "seeing" means 
that if every point of s ; sent out horizontal light rays then some 
of them would hit Sj. (Note that the relation is symmetric, even 
though this description isn't). In the example collection shown 
in figure 26, 52 and s 4 see each other, while si and do not. 
This visibility relation between segments can be useful to us 
because, if we remove from S a collection of mutually invisible 
segments to obtain a reduced set S', then the answer to the left- 
and-right neighbor problem in S' is not too different from that 
in S. Specifically, if we think of the deleted segments as having 
become transparent, then the rays we emit from our query point 
can only pass through at most one transparent segment before 
they encounter some segment of S' (or go off to infinity). 

Assume that we have preprocessed our collection of seg- 
ments as mentioned above so we have for each segment endpoint 
its left and right neighbors. It is most useful to record this in- 
formation on the neighbors: for each segment s we maintain 
two lists ordered in y that describe the segment endpoints for 
which s is the left or right visible segment, respectively. These 
are called the visibility lists for s and their total length is called 
the degree of s. We require that a left-and- right neighbor query 
from a point p report not only the two segments stopping the 
rays from p but also the specific interval in the right visibility 
list of the left neighbor and the left visibility list of the right 
neighbor where these rays fall. Note that the total size of the 
visibility lists, summed over intervals, is at most 4n; further- 
more, these lists can be set up in 0(n log n) time by sweeping 
S horizontally. 

Now the following sampling process suggests itself: remove 
a large collection of mutually invisible segments from S to ob- 
tain a coarsened subcollection S'. For this subcollection S' we 
store separately the original visibility lists in S of the segments 
in S'. We annotate these lists by giving for each segment s in 
5', and for each pair of adjacent y values in the left and right 
visibility lists of s, the name of the neighboring deleted seg- 
ment, if any. This is the deleted segment that lies immediately 
to the right or left of the appropriate interval in the lists for s. 
We refer to these auxiliary structures as the support of S'. If we 
have the support of S', then we can take the answer to a left- 
and-right neighbor query in S' and compute the corresponding 
answer in S in constant extra time: we need to check only the 
single deleted segment of 5 that might have gotten in the way 
in each direction. See figure 27. 



S' 



Figure 27. 
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From the visibility lists for 5 we can derive the visibility 
lists for S' by a pretty straightforward process. When segment 
s is to be deleted, its left visibility list has to be merged with 
those of its right neighbors, and correspondingly for the right 
visibility list. These neighbors are visible from s, so they are not 
themselves deleted. The merging process can be carried out in 
time proportional to the sum of the lengths of the two visibility 
lists for s. Once we have the proper visibility lists for S' itself we 
are back where we started, so we can repeat the whole sampling 
process over again. 

Selecting the sample. For this approach to work we have to 
make sure that we can always find a sufficiently large set of 
mutually invisible segments to remove, and then do the list 
updating efficiently. We have already remarked that the total 
size of all the visibility lists for all segments is less that 4n. Thus 
the average segment has degree less than 4. Equivalently, there 
are at least n/2 segments of degree less than 8; let us call their 
set T. From among the segments in T we would like to extract a 
large set of mutually invisible segments. So let us take the first 
segment from T, place it in a bag, and then delete from T each 
of the at most 7 other segments that it sees. Now continue this 
process till no more segments are left in T. Since at each step 
at most 8 segments are removed altogether, there will certainly 
be at least n/(2 X 8) = n/16 mutually invisible segments placed 
in the bag. These segments form the set that we remove from 5 
to obtain 5'. Note that all the removed segments have degree 
less than 8. 

The process of computing the segments to be removed from 
S clearly takes linear time. The auxiliary support structure for 
S' can also obviously be stored in linear space. Now the set 
S' is at most 15/16 of the size of 5. Furthermore, its visibility 
structures can be computed from those of S in time linear in n, 
since all the merging steps that have to be done involve constant 
size lists. 

The hierarchical structure. Now our solution is clear: we 
iterate this sampling process till we get to some constant-size 
collection of segments R. We store the support structures for 
all the intermediate samples up to and including R. The cost of 
computing and storing each level is proportional to the number 
of segments in it, which is at most a constant fraction p = 15/16 
of the segments in the previous level. Therefore, the total cost 

(time and space) of all these samples is 0(n + pn + p 2 n H ) = 

0(n/(l-p)) = 0(n). 

In order to answer a left-and-right neighbor query we use 
an exhaustive search on R. Then, using the support structures 
we have saved, we back up through this hierarchy of samples, 



Geometric Algorithms 



29 



at constant cost per level. From the previous discussion it is 
apparent that the number of levels is O(logn), which is also 
the cost of answering a query. We have shown the following: 

Theorem 13. It is possible to preprocess n vertical segments in 
0(n log n) time and 0(n) space so that any left-and-right 
neighbor query can be answered in 0(log n) time. 

Rectangular point location. A minor variation of the tech- 
nique presented here can be used for solving the problem of 
point location in a rectangular subdivision of the plane. This 
problem has obvious applications to VLSI design tools and 
multi-window user interfaces. We look at the vertical sides of 
our rectangles and then lump them together into maximal ver- 
tical segments. The rectangular partition denned by these ver- 
tical segments is almost the original subdivision. It only fails to 
account properly for vertical towers of rectangles with identical 
width. However, those are not difficult to incorporate into the 
technique shown. The result is a point location algorithm as 
asymptotically efficient as that of Kirkpatrick, but which stays 
entirely in the rectangular domain and does not require any 
coordinate operations other than simple comparisons. 
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Part B: Data Structures 

In this part we focus on a number of data-structuring tech- 
niques that have been found especially useful in designing effi- 
cient geometric algorithms. 

Bl. Fractional cascading 

In computational geometry many search problems and range 
queries can be solved by performing an iterative search for the 
same key in separate ordered lists. In this section we demon- 
strate a technique known as fractional cascading [10] for effi- 
ciently solving such problems. Fractional cascading is a data- 
structuring technique for solving the iterative search problem, 
which we formulate as follows: let G be a directed graph whose 
vertices are in one-to-one correspondence with a set of sorted 
lists; given a query consisting of a key q and a connected sub- 
graph 7T of G, search for q in each of the lists associated with the 
vertices of 7r. This problem has a trivial solution involving re- 
peated binary searches, at a logarithmic cost per list. Fractional 
cascading establishes that it is possible to do much better: un- 
der the assumption that the graph G has bounded degree, we 
are able to organize the set of lists so that we have to pay the 
logarithmic cost only once, when searching the first list, and 
then pay only a fixed amount to look up q in each of the re- 
maining lists. This can be done using only a linear amount of 
extra storage. These bounds are clearly best possible. 

The interval stabbing problem. We now consider the fol- 
lowing problem: given a collection of n intervals on the real 
line, organize them into a data structure so that the number 
of them containing an arbitrary query point can be efficiently 
computed. We would like a method that uses linear space and 
solves this problem in time O(logn). In the sequel we assume 
that the given intervals have distinct endpoints. 

We organize our intervals into a tree T using an idea due 
to McCreight [45,46]. First we normalize all the endpoints so 
that they become the rationals l/(2n + 1) to 2n/(2n + 1), in 
left-to-right order. These rationals are all in the interval (0, 1). 
We put in the root of T the collection of intervals containing the 
rational 1/2. In the left child of the root we put the intervals 
not already placed that contain the rational 1/4, and in the 
right child of the root we put those not already placed that 
contain 3/4. Note that any interval containing both 1/4 and 
3/4 already contains 1/2, so it has been allocated in the root 
node. We continue in this process until every interval has been 
allocated to exactly one node of the tree T thus formed. The 



32 



Geometric Algorithms 



normalization assumption implies that the depth of T will be 
only O(logra). 

Each node v of T thus ends up containing a collection of 
intervals. We construct two ordered lists from these intervals 
and store them in v. One list contains all the left endpoints 
of the intervals in v sorted by a; value, and the other the right 
endpoints, similarly sorted. Clearly the total space that the tree 
structure takes, together with these ordered lists, is 0(n). 

Processing a query. A query asking for the number of inter- 
vals containing some point p can then be answered as follows. 
We start at the root and traverse a path down the tree, essen- 
tially emulating a binary tree search for p, where the key value 
of the root is 1/2, of its left child 1/4, etc. For each node v 
visited we compare p to its key value. If p is greater than the 
key value, then we locate p by a binary search among the right 
endpoints of the intervals stored with v. Once we know the lo- 
cation of p, we can easily obtain the number of intervals in v 
containing p. Then we continue the search with the right child 
of v. We perform the symmetric operations if p turns out to be 
less than the key value in v. This method works because once 
we have (say) branched right from a node v we can be sure that 
no interval in a node that is in the left subtree of v can possibly 
contain p. 

The cost of this method is 0(log 2 n): we have 0(log n) levels 
to traverse, and at each level a binary search to perform. How 
could we reduce the O(log 2 n) overhead to only O(logn)? 

Using fractional cascading. Here is where fractional cascad- 
ing comes in. The underlying graph G is the tree T. The sorted 
lists are the ordered end point lists (we treat the list of left end- 
points entirely separately from the list of right endpoints). The 
degree is bounded by three. So here is what we can do. Consider 
first the left endpoint lists only. Each leaf of T extracts every 
other one of its list elements and merges them into the list of 
its father. A node with two leaves as children has its own list 
merged with that of the "sampled" lists from each of its two 
children. Now each of these fathers with two children leaves in 
turn samples his own new list by taking every other element 
and sends that sample up to be merged with the list of his own 
father. This process continues all the way up the tree: whenever 
a node has merged with samples from both of his children, he 
samples himself and sends the sample to his father. Because 
of the acyclic structure of the tree T this process eventually 
terminates. When it has done so, each node of v contains an 
augmented list of left endpoints. This list for v still contains all 
the left endpoints of intervals stored in v, but it also contains a 
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sample of those in its children, a smaller sample of those in its 
grandchildren, and so on. 

Analysis. First of all it is clear that if we know where we are 
in the augmented list of a node v, then we can easily compute 
where we are in the regular list of v. Second, if we know where 
we are in the augmented list for v, we nearly know where we 
are in the augmented lists of either of its children: at most 
one comparison with an unsampled value in the interval that 
we are in will tell us our position exactly in a child. This is 
an important accomplishment, because once we have done a 
binary search to locate p in the augmented list of the root, we 
can then descend the tree T and locate p in the relevant lists 
at only constant additional cost per node visited. Thus this 
propagation of samples has made neighboring lists sufficiently 
similar that the search overhead is reduced to only 0(log n): an 
initial binary search, and then constant cost per level. 

But what about storage? Hasn't all the propagation of sam- 
ples increased storage use significantly? The answer is no. For 
each list, half of it went to its father, half of that went to the 
grandfather, and so on. So that the total size of all the samples 
of a particular list floating around is a geometric series in the 
original length of the list, and therefore linear in that length. 
This argument is not quite rigorous, but it conveys the intuition 
why this sampling process has in fact not more than doubled 
the storage used. For a rigorous argument the reader is referred 
to the paper by Chazelle and Guibas [10]. This notion of prop- 
agating cascaded fractions of the whole forms the essence of 
the idea behind fractional cascading. Thus we have shown the 
following: 

Theorem 14. In 0(n log n) time it is possible to organize n 
intervals into a linear space data structure such that the 
number of intervals containing a query point can be found 
in time O(logn). 



B2. Finger trees 

A classical data structure that has found a number of im- 
portant uses in computational geometry is the finger search tree. 
A finger search tree is like an ordinary balanced tree used to 
implement the dictionary operations of find, insert, and delete. 
Some additional structure is provided, however, so that these 
operations are especially efficient in the vicinity of certain pre- 
ferred positions, indicated by fingers. Let n be the size of the 
dictionary; while in an ordinary balanced tree each operation 
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Figure 28. 



requires time O(logn), in a finger search tree the time is only 
O(logd), where d is the distance from the finger where the 
search begins. Thus if there is locality of reference, by keep- 
ing a finger in the active region of the key space we can speed 
up all three operations. The finger search tree structure was 
introduced by Guibas, McCreight, Plass, and Roberts [30] and 
has since been further developed by many researchers. 

For our discussion here we will need the version described 
by Tarjan and Van Wyk under the name of homogeneous finger 
search tree [72]. The structure consists of a 2,4-tree (or, equiva- 
lent^, a red-black tree) with additional internal links from each 
node to its parent, and to its left and right neighbors on the 
same level of the tree. See figure 28. These extra pointers 
allow fast searches starting from any given node. More specif- 
ically, the cost of locating an item that is d items away from 
a given node (in the linear order of nodes) is only O(logtf). In 
fact, if we maintain two fingers pointing at the first and last 
leaves of the tree (or, alternatively, if we make the left /right 
neighbor links circular in each level), we can reduce this to 
0(logmin{d, n — d}), where n is the total number of nodes in 
the tree. To achieve this bound, we start searching simultane- 
ously in both directions from the given node, wrapping around 
when we reach either end, and stopping both searches as soon 
as one of them succeeds. 

Homogeneous finger search trees also allow us to perform 
efficiently a few other list operations, such as insertion, dele- 
tion, splitting, and concatenation. In particular, once we have 
located the position of a given key in the tree, we can insert 
or delete the corresponding node in only 0(1) additional time. 
Also, given any pair x, y of leaves, we can split the tree into two 
smaller trees, with all the nodes between x and y in one tree and 
the balance in the other, in time 0(log min {d, n - d}), where 
d is the number of nodes between the two leaves. Conversely, 
if we are given two trees of sizes d and n - d whose keys lie 
in disjoint segments of the key space, we can join them into a 
single tree in the same time bound. For a careful description of 
these algorithms see the appendix to the triangulation paper of 
Tarjan and van Wyk [72]. 

Note that these are amortized bounds. What this means 
is the following. Consider any sequence of searches, insertions, 
deletions, splits and joins starting with an empty tree. Let n t - 
be the total size of the operands for the ith operation, and d t be 
either the distance between the starting node and the affected 
node (for find, insert, delete), or the size of the first operand (for 
join), or the size of the first result (for split). What the amor- 
tized bounds say is that the total cost of the whole sequence 
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will be O(£\logmin{d,-,ni - d,}), even though the ith opera- 
tion by itself may cost much more than O(logmin {di, n,i - e?<}). 
We will discuss amortized analysis in more detail in section C3. 

Jordan sorting. The application which we will use to demon- 
strate the use of finger search trees is known as Jordan sorting. 
Suppose we are given a simple polygon P of n sides and an 
(infinite) straight line /. We wish to compute the intersections 
of P and /, sorted in the order in which they occur along /. It 
is possible to do this in only 0(n) time, by an elegant algo- 
rithm due to Hoffman, Mehlhorn, Rosenstiehl, and Tarjan [36]. 
In what follows we outline their solution. 

Consider what happens to one side of /; the other side is 
symmetric. The line / breaks up the polygon into a bunch or 
arcs. See figure 29. These arcs are either disjoint or nested, 
since the polygon is simple. We represent the containment re- 
lationships by a tree S, as shown. Each node is an arc and its 
children are all arcs directly nested within it, sorted by the order 
in which they occur along /. The plan is to follow the polygon 
P, and, as each arc of P on the current side of / is completed, 
update the tree S. We do this on both sides of / simultaneously, 
with pointers connecting every arc endpoint on one side to the 
corresponding arc endpoint on the other side. When we are fin- 
ished we can just make a symmetric order traversal of S and 
derive the sorted order of intersections of P with /. 

So how do we update 5? Whenever we start processing a 
new arc 7, we look for the innermost pair of previously processed 
arcs a, a' whose endpoints bracket the starting point of 7. As we 
will see, while processing the previous arc (on the other side of /) 
we will have determined the nearest arcs that bracketed its exit 
point; therefore, we can get a and a' in 0(1) time by following 
the pointers mentioned above. If the neighbors are nested, that 
is, they are parent and child in 5, then the outermost of the 
two will become the father ip of the new arc 7. If they are not 
nested, they must be siblings in 5, and their common father <p 
will become the father of 7. See figure 30. 

Next, we follow 7 until it reaches / again. The point where 
this happens must lie in the gap between two children /?,/?' of 
<p, or between tp and its first or last child. We locate (3 and 
/?' in the children list, and then we update S by inserting the 
new node 7 as a child of <p, and demoting all children of (p that 
are enclosed by 7 (a through /3 in the figure) to be children of 
7. In order to do these operations efficiently, we represent the 
children list for each node of S as a homogeneous finger search 
tree. (These tree structures should not be confused with the 
tree S. Rather, we should think of each finger tree as a linear 
list of children, with some hidden machinery that allows us to 
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efficiently perform certain search and update operations on the 
list.) 

Analysis. It remains to analyze the complexity of this algo- 
rithm. Observe that the insertion of each arc 7 requires locat- 
ing its exit point among the children of <p, splitting off all the 
enclosed children (if any), and inserting 7 in their place. The 
problem is finding an upper bound to the total cost of these 
operations. 

Again, let's consider only what happens on one side of the 
line /, and let S be the final nesting tree for this side. Now con- 
sider the situation after the algorithm has just added some arc 
(p to the tree. It is helpful to think that all the arcs in S are 
initially "penciled in" in the plane, and then redrawn in ink in 
the order in which the algorithm inserts them. Suppose cp has m 
descendants, excluding itself, in the final tree S. In other words, 
<p encloses m other arcs, penciled or inked, arbitrarily nested. 
Suppose also that of these descendants A; are still penciled in. 
We call the quantity c = m + k the complexity of <p. We wish to 
determine the maximum cost T(c) (as a function of the com- 
plexity c) for completing the subtree rooted at <p, that is, for 
adding the arcs that are descendants of <p in 5 but haven't been 
processed yet. Note that some of the descendants of <p may have 
been inked before (p; we aren't interested in how much their in- 
sertion cost, nor in how much it cost to insert <p — we just care 
about the cost of completing the inking process. Note also that 
the cost of inking the arcs under <p is not affected by any arcs, 
inked or not, that lie outside <p; therefore, we can ignore them 
completely in the determination of T(c). 

Now let 7 be the first descendant of (p to be added after <p; 
suppose that 7 has p, p < m, descendants, of which q, q < k, 
are still penciled in. The complexity of 7 is then d = p + q < c. 
Since we are representing the inked children of ip by finger trees, 
if (p currently has i inked children, and 7 covers j of them, 
the (amortized) cost of adding 7 is O(log min { j + 1, i - j + 1}). 
Obviously, the number j is no greater than the total number p 
of descendants of 7 in 5. Similarly, i - j + 1 is no greater than 
m-p, the descendants of <p that are not descendents of 7. Since 
P < P + 9 = d, and m - p < (m + k) - (p + q) = c - d, the cost 
of adding 7 is at most O(logmin {d + l,c - d}). 

Once 7 is inked, what happens to the arcs below 7 is inde- 
pendent of what happens to those between 7 and tp. The cost 
of finishing the former is by definition at most T(d). The cost 
of finishing the latter would be the same if the descendants of 
7 did not exist, i.e. it is as if <p enclosed only m-p arcs, of 
which k-q-1 were still penciled in (the -1 is there because 
we just inked 7). Since (m-p) + (k-q- 1) = c - <f - 1, this cost 
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is at most T(c — d — 1). Therefore we can write the following 
recurrence for the cost T(c): 

T(c) < max (T(d) + T(c - d - 1) + 0(logmin {d + l,c - d})). 

The solution to this recurrence is T(c) = 0(c) [49, vol. I]. Since 
the complexity of the first arc is less than twice the number of 
arcs n, the cost of building the entire tree S is only 0(n). We 
conclude the following: 

Theorem 15. The intersections of a simple polygon and a straight 
line, in the order in which they occur along the line, can be 
computed in linear time. 



B3. Persistent data structures 

When we discussed the plane sweep technique in section A4, 
we observed that it transforms a static two-dimensional problem 
into a dynamic one-dimensional one. Objects enter the active 
list as the sweep line hits them, and exit the list as the sweep 
line moves past them. In some sense, the history of the active 
list represents all the data present in the original problem. Some 
geometric problems can be solved efficiently by encoding this 
history in a compact data structure that allows fast access to 
any of the active lists recorded during the sweep. Noe that this 
is in a sense the opposite of the plane sweep paradigm, for we are 
transforming a dynamic one- dimensional problem into a static 
two-dimensional structure. 

In general, a persistent data structure is one that accepts an 
arbitrarily long sequence of updates, but is able to remember 
at any time all its earlier versions. Driscoll, Sarnak, Sleator, 
and Tarjan have given general techniques for taking ordinary 
ephemeral data structures and making them persistent [19]. An 
especially useful example is that of persistent sorted sets. In 
such a set we maintain a linearly ordered set of items. Items are 
always inserted and deleted from the "last version" of our sorted 
set. Queries as to the position of a particular item, however, can 
be made either in the present or in the past: we are allowed to 
ask for the place of an item in the sorted list as it was in some 
earlier version. If the current list contains n items and a total 
of m updates have been made, then a structure proposed by 
Sarnak and Tarjan [62] achieves update time of O(logn), access 
time to any version in the present or past in 0(log m) time, and 
needs 0(1) amortized space per update, starting from an empty 
set. 
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Application to point location. Persistent sorted lists can be 
the basis of an alternative point-location algorithm for planar 
subdivisions. Let 5 be a subdivision of the plane into polygonal 
regions with a total of n edges. We assume that no edges of 
S are vertical. Let's cut the plane into a sequence of vertical 
slabs, by drawing a vertical line through every vertex of S. 
Within each slab 5 looks like a stack of trapezoids and triangles 
separated by line segmens that cut completely through the slab, 
and thus have a well-defined top-to-bottom order. See figure 31. 
If we could afford to represent the segments in each slab by a 
separate sorted list, then we could locate a given query point 
x in O(logn) time, by first locating the slab containing x, and 
then locating the two consecutive segments in that slab that 
bracket the point x. 

Unfortunately storing the slabs separately would require 
0(n 2 ) storage in the worst case. However, we can reduce the 
space by noting that the sorted sets of line segments that form 
the trapezoid bounding edges in contiguous slabs are similar. 
Consider sweeping the plane from left to right by a vertical 
line. While the sweep line is inside a slab, the active list contains 
precisely the segments that cross the slab. As we sweep over the 
boundary from one slab to the next, some segments are deleted 
from the activel list and some are added. Over the entire sweep, 
there will be 0{n) insertions and deletions, one insertion and 
one deletion per edge of 5. Therefore, if we represent the active 
list as a persistent sorted set of segments, by the end of the 
sweep this data structure will use only O(n) space, and will 
allow us to search any slab in O(logn) time. 

Theorem 16. Given a polygonal subdivision S with n edges, we 
can build in 0(n log n) time and 0(n) space a search data 
structure that allows us to find the region of S containing 
an arbitrary query point x in only O(logn) time. 
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Part C: Analysis Techniques 

For the major part, the techniques used for the analysis of 
geometric algorithms are not intrinsically different from those 
used in the analysis of algorithms in general. It is often the 
case, however, that some of the quantities needed in the analy- 
sis depend on non-trivial results of topology and combinatorial 
geometry. Below we present two such results, having to do with 
the length of Davenport-Schinzel sequences (section C2), and 
the total complexity of the faces incident to a line in a line ar- 
rangement (section CI), which have been used recently in the 
analysis of many geometric algorithms. 

In the theoretical kind of computational geometry that we 
are considering in this paper, the efficiency of an algorithm is 
most often measured by its asymptotic worst-case time and 
space requirements, viewed as a function of the input size, out- 
put size, or some similar measure of the problem's complexity. 
In fact, the large variability in the output size of many geomet- 
ric algorithms makes this quantity an important parameter of 
the analysis. The reader is referred to Seidel's dissertation [64] 
for many examples of output-sensitive algorithms. An increas- 
ingly common tool is amortized analysis, which we already used 
in section B2 to derive a linear bound for Jordan sorting. In sec- 
tion C3 we show a much simpler application of this technique 
to the analysis Graham's convex hull algorithm. 

Average-case analysis of geometric algorithms is rarely done, 
and few general-purpose techniques have been developed for 
this area. Section C4 shows one of the few results obtained so 
far, an asymptotic analysis of the expect number of vertices on 
the convex hull of n points uniformly distributed in a triangle. 

CI. Exploiting results from combinatorial geometry 

In the analysis of algorithms in computational geometry 
one often must make use of results of a combinatorial nature 
dealing with geometrical entities. Such results properly form the 
realm of a field known as combinatorial geometry. Griinbaum's 
book on convex polyhedra [28] is a good example of a work in 
this area. 

We illustrate the use of such techniques in constructing the 
arrangement of n lines in the plane, that is, the subdivision of 
the plane determined by n given lines. An interesting combina- 
torial fact about such arrangements is the following: suppose we 
look at all faces of the arrangement adjacent to one of the lines 
in the arrangement, the so-called horizon of the line. Then the 
total number of edges on all these faces is only O(n). Notice 
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Figure 32. 




Figure 33. 



that even one face can have n edges, so this bound is not as 
trivial as it may seem. 

Based on this combinatorial fact we can design an algorithm 
for building the arrangement of the n lines by adding them in 
one at a time. In order to add a line /, we start at the appropriate 
infinite face R (one of the two, actually) and find the unique 
edge e through which / exits R. Then we cross over to the region 
on the other side of e and proceed in the same way. At each step, 
we find the edge where / exits the current region by simply 
walking along the boundary of the latter and testing each edge 
against /. Clearly, the total cost of this process is proportional 
to the total number of edges on the line's horizon, which as we 
have said above is only 0(n). Therefore the total cost (in both 
time and space) for inserting all n lines is 0(n 2 ). If our aim 
is to build and save the arrangement, then this is clearly best 
possible, as the arrangement will contain 0(ra 2 ) vertices, edges, 
and regions if the lines are in general position. This algorithm 
was given by Chazelle, Guibas and Lee [12], and Edelsbrunner, 
O'Rourke and Seidel [22]. 

Now we prove the combinatorial result mentioned above. 
We assume that the arrangement contains no parallel lines. Let 
/ be the chosen line. Consider first the faces of the arrangement 
that are adjacent to and above /. Since the other n — 1 lines cut 
/ into n segments, there are exactly n such faces. For each face 
/, let v be the vertex furthest from Z; That point is unique, by 
our non-parallelism assumption. Call the edge that lies on I the 
floor of /. The floor and the vertex v separate the remaining 
edges into the left and right sides of /. The edges adjacent to 
v are the ceiling of /. (If the face is unbounded, its ceiling is 
by definition the one or two edges that go off to infinity). The 
edges that are neither floor nor ceiling, if any, are called the 
walls of /. See figure 32. 

We claim that each line of the arrangement can occur in 
the wall of at most two faces, once as left wall and once as right 
wall. To see why, suppose a line to occurred twice as a right 
wall, on two faces / and g. Without loss of generality, suppose 
/ is the one closer to the intersection of m and I. See figure 33. 
Let to' be the next line on the boundary of /, counterclockwise 
from to (which may be either a wall or a ceiling of /). Since to 
and to' are on the right side of /, we conclude m' meets I to the 
right of m. Hence, to the left of its intersection with m the line 
m' lies always between I and m. Therefore it must cut through 
the face g, a contradiction. 

A similar argument holds for left wall occurrences. We con- 
clude that the total number of wall edges in all regions under 
consideration is at most twice the number of lines, that is, 2n-2. 
Actually, the total is at most 2n - 4, since the line that crosses 
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/ furthest to the left cannot appear as a left wall, and symmet- 
rically for the right. Moreover, each face has one floor edge and 
at most two ceiling edges (with the first and last regions hav- 
ing only one). The total number of edges in all those regions is 
therefore at most 5n - 6 = 0(n). (This result can be proved 
also by a Davenport-Schinzel argument, or through the use of 
Theorem 17; details are left to the reader.) We conclude that 

Theorem 17. The arrangement of n lines in the plane can be 
built in 0(n 2 ) time and space. 



C2. Davenport-Schinzel sequences 

Davenport-Schinzel sequences are interesting and power- 
ful combinatorial structures that arise in the analysis and cal- 
culation of the lower envelope of collections of functions, and 
therefore have applications in many geometric problems that 
can be reduced to the calculation of such an envelope. Roughly 
speaking, an (n,s) Davenport-Schinzel sequence is a sequence 
composed of n symbols with the properties that no two adjacent 
elements are equal and that it does not contain as a subsequence 
any alternation of two distinct symbols of length 3 + 2 (e.g. an 
(n,3) sequence is not allowed to contain any subsequence of 
the form a ... b ... a ... b ... a). The main goal in the analysis of 
these sequences is to estimate their maximal possible length for 
any given values of the parameters n and s. 

The importance of Davenport-Schinzel sequences lies in 
their relationship to the combinatorial structure of the lower 
envelope of a set of functions. Let /i,/2,...,/n ben real- valued 
continuous functions defined on the real line, having the prop- 
erty that each pair of them intersect in at most s points. The 
graph of their lower envelope /, defined by f(x) = min t - /,(x), 
is the concatenation of a finite number of arcs, where the kth 
arc belonging to some function fi h . The sequence of indices 
ii,i - 2,...,ifc from left to right is an (n,s) Davenport-Schinzel 
sequence. Conversely, any (n, s) Davenport-Schinzel sequence 
can be realized in this way for an appropriate collection of n 
continuous univariate functions each pair of which intersect in 
at most s points. 

The crucial property of these sequences is that, for a fixed 
s, the maximum length A s (n) of an (n,s) sequence is "prac- 
tically linear" in n. More precisely, Ai(n) = n, A2(n) — In — 
1, and \j,{n) = 0(na(n)), where a(n) is an extremely slow- 
growing function of n, the functional inverse of Ackermann's 
function. Similar bounds have been proved for A a (n) when s > 
3 [34,67,69]. 
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This property makes Davenport- Schinzel sequences a use- 
ful and powerful tool for solving numerous problems in dis- 
crete and computational geometry. Indeed, the original moti- 
vation for the introduction of these sequences by Davenport 
and Schinzel was a geometric problem arising in analysis of the 
combinatorial properties of the pointwise maximum of a collec- 
tion of solutions to a linear differential equation. In the last few 
years these sequences have been extensively studied by Micha 
Sharir and various coworkers. They have found application to 
a diverse set of geometric problems, including an analysis of 
the exterior boundary of the union of (intersecting) segments 
in the plane, the calculation of Euclidean shortest paths in 3- 
space amidst polyhedral obstacles, the estimation and compu- 
tation of the configuration space used in solving various motion 
planning problems amidst polygonal obstacles, the analysis of 
certain time-varying geometric configurations, and many oth- 
ers [3,63,68,70]. 

Complexity of red-blue intersections. We will illustrate the 
use of Davenport-Schinzel sequences on the following problem. 
Suppose we are given two sets of convex polygons, "red" and 
"blue," such that no two polygons with the same color intersect. 
See figure 34. A red and blue pair may have points in common, 
in which case their intersection is a "purple" convex polygon. 
Observe that the purple polygons are pairwise disjoint. Let the 
red and blue polygons have total size m and n, respectively 
(where size = number of edges). Now suppose we take a fixed 
number k of purple polygons. Our problem is to estimate the 
maximum total number of edges <r(m, n, A;) of these k polygons. 

The difficulty here is that a red or blue edge may give rise 
to many purple edges. It is easy to see that a purple region P 
can have m + n sides (exactly) in the worst case, since each red 
or blue edge contributes at most one edge to P. On the other 
hand, the number of purple polygons can be as high as 0(mn) 
in the worst case (consider for example a grid of m/4 vertical red 
stripes and n/4 horizontal blue ones), but then their total size 
is only 0(mn) instead of 0(mn(m + n)). That is, one purple 
region may have many edges, but if we take a lot of regions 
their average size will be quite small. For intermediate values 
of k, the answer is far from obvious. For example, figure 35 
shows that ct(6,3,2) > 11. Can we get a(m,n,2) > c(m + n) 
for some c > 1? Surprisingly, the answer is no. It turns out that 
<r(m, n,k) < m + n + 4(k - 1) for all m, n > 3 and k > 0. 

We will now prove this result. To avoid confusion, let's call 
the purple regions cells and their sides arcs, and reserve the 
name edge to those of the red and blue polygons. So, let K be 
some set of k cells. Let us fix our attention on some particular 
red polygon Pi with m,- edges. Let K{ be the set of those cells 
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from K that are contained in Pi, and ki the number of such 
cells. See figure 36. Let us classify the arcs of Ki as "red" 
or "blue" depending on their origin. We wish to show that the 
number of red arcs in K{ is at most m f + 2(Jb,- - 1). 

Label each red arc of K t with the name of the unique blue 
polygon that contains it. Now consider the sequence of labels we 
encounter as we go once around the polygon Pi. Of course sev- 
eral arcs may carry the same label, since two convex polygons 
may weave in and out of each other arbitrarily often. However, 
between any two consecutive arcs with the same label the poly- 
gon Pi must have at least one vertex. We "charge" the second 
of the two arcs to that vertex, and delete the corresponding 
entry from the label sequence. Clearly, every vertex of Pi will 
be charged at most once, which means at most m t - arcs are 
accounted for this way. 

How many red arcs may still be unaccounted for? Here is 
where the Davenport-Schinzel theory comes to the rescue. Con- 
sider the sequence of labels remaining after the deletions above. 
Each letter may still appear several times, but not twice in a 
row. Moreover, we cannot have two pairs of repeated letters in 
the pattern ...a. a. ..6.... If that occurred, then convex- 
ity would imply that the blue polygons a and 6 intersect, which 
is a contradiction. See figure 37. This means the string of un- 
accounted red edges is a circular Davenport-Schinzel sequence 
on ki letters with parameter s = 2. As mentioned earlier, the 
length of such a sequence is at most 2(fc t - - 1). We conclude Ki 
has a total of at most m< + 2(ki — 1) red arcs. 

If we sum this bound over all red polygons Pi that contain 
at least one cell we conclude that the number of red arcs in the 
k given cells is at most m + 2(k - 1). Similarly, the number of 
blue arcs is at most n + 2(k — 1). We have therefore proved the 
following claim: 

Theorem 18. The total size of any k of the cells formed by 
intersecting two collections of disjoint convex polygons of 
total size m and n is at most m + n + 4(k — 1). 




Figure 36. 




C3. Amortized analysis 

We already encountered the concept of amortized complex- 
ity in section B2, in the analysis of finger search trees. The tech- 
niques of amortized analysis [71] have also found a number of 
interesting uses in the analysis of geometric algorithms. Perhaps 
the earliest significant algorithm in the field, Graham's method 
of computing the convex hull of n points in the plane [26] al- 
ready offers an example of amortized analysis. In that method 
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one first sorts the n points polarly axound some point that is 
interior to the convex hull. Following that, one performs what 
has come to be called the Graham scan. This is a pass over the 
points in polar order where we maintain a convex chain as we 
go along; this chain is our best guess as to what the convex 
hull is, based on points we have seen so far. If the next point 
considered destroys the convexity of the chain, then one has to 
back up through the previous points in the chain and discard 
points till convexity is restored. This is most conveniently done 
by keeping the current convex chain as a collection of points in 
a stack. If the new point fits the chain, then it is just stacked 
as well. If, on the other hand, the new point destroys the con- 
vexity, then points are successively popped from the stack till 
convexity is restored between the stack contents and the cur- 
rent point. That point is then pushed onto the stack and the 
process repeats. Since each point is stacked and unstacked at 
most once, the whole process takes time linear in n. 

Notice that the processing time for each point during the 
Graham scan can be highly variable. Some points may use 0(1) 
time, and others may take Q(n) time. Still, the charging scheme 
of putting the cost of the operations on the points themselves as 
they are stacked or unstacked and not on the points that cause 
the stacking or unstackingis an example of amortized complex- 
ity analysis. Such arguments frequently involve clever charging 
schemes, such as the use of potential functions. Another exam- 
ple where amortized complexity comes in is the analysis of the 
topological plane sweep method for an arrangement of n lines 
by Edelsbrunner and Guibas [21]. 

C4. Average-case analysis 

The field of average-case analysis of geometric algorithms 
is less well developed than its counterpart in combinatorial al- 
gorithms. One difficulty is that most geometric problems have 
infinite input domains, for which there are no natural probabil- 
ity distributions. For example, what does it mean to choose n 
random points on the plane, or a random simple polygon with 
n vertices? Moreover, any given input distribution is unlikely to 
be represetative of data to be encountered in practice: the sim- 
ple polygons typically encountered in cartography applications 
are statistically quite different from those typical of engineering 
design. Finally, the non-linearity of most geometric operations 
causes internal results to have extremely complicated distribu- 
tions. Just to mention a simple example, suppose we take two 
points p and q with random coordinates, independently drawn 
from a Gaussian probability distribution. The coefficients of the 
line joining p and q will be second-degree rational functions of 
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all four coordinates; as such, they will not be independent of 
each other, and they will not follow a Gaussian distribution, or 
any of the distributions normally studied in statistics books. 

The study of this kind of problems is known as integral or 
stochastic geometry. For an introduction to this field, see the 
work of Santalo [60,61] and Miles [50,51,52,53]. We expect 
that techniques from this field will find more applications in 
the analysis of geometric algorithms in the future. 

Random convex hull. As an example of average-case anal- 
ysis, we will derive an asymptotic expression for h(n), the ex- 
pected number of vertices of the convex hull of n random points 
in the plane. The importance of this quantity stems from the 
fact that many geometric algorithms start by computing the 
convex hull of n given points, and often the cost of subsequent 
steps is determined by the size h of the hull, rather than the 
total number of input points. 

The value of h(n) obviously depends on the probability dis- 
tribution of the given points, as well as on their number. In this 
section we will consider a very simple case: we will assume the n 
given points are uniformly and independently distributed inside 
a triangular region of the plane. Following the approach used 
in 1963 by Renyi and Sulanke [59], we will prove the following 
result: 

Theorem 19. The average number of vertices in the convex hull 
ofn random points uniformly and independently distributed 
in a triangle is 



Proof: First, observe that affine maps preserve convexity, "in- 
side-outside" relations, and uniform distributions. Moreover, for 
any two triangles there is an affine map that takes one into the 
other. We conclude that if the theorem holds for the equilateral 
triangle T with unit side, it will automatically hold for any 
other triangle. 

Let P = {pi,P2, • ■ ■ ,Pn} be the input points. If we define 
the random variables 



h(n) = 2H n - X 



= 2(1+ _ + _ + . .. + __) 



= 2(lnn + 7) + o(l) 



(1) 




1 if piPj is an edge of hull(P), 
0 otherwise. 



(2) 
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then the number of edges (vertices) of hull(P) will be 

l<i<j<n 



(3) 



and therefore 

h(n)= £ Pr ( £ u = 1 ) ( 4 ) 

l<«<i<" 

Since all points are drawn from the same distribution, the prob- 
abilities Pr(eij=l) are all equal, and it suffices to compute one 
of them, say Pr(£i 2 = 1). This number is the probability that 
points P3, p4, . . . p n fall on the same side of the line P1P2 . 

With probability one, the line pip 2 divides the triangle T 
into a smaller triangle T\ 2 and a convex quadrilateral Q\ 2 . See 
figure 38. Let u and v be the lengths of the two sides of T\ 2 
that are not on the line pip 2 . The areas of T and T\ 2 are then 



1 . 7T 

m = -„„- = - 

. . 1 . 7T >/3 

Xi2 = —uv sin — = uv — 
1 1 2 3 4 



(5) 



Therefore, |Ti 2 | / \T\ = uv, and \Q n \ /\T\ = l- uv. The prob- 
ability Pr(£i2=l) of all remaining points falling on one side of 
P\P2 is then 



J J {{uv) n - 2 + {l-uv) n - 2 )Vr{p l )dp 1 Yt{p 2 )dp 2 



(6) 



T T 



where Pr(p) is the probability density of the input points at p, 
namely 1/|T| = 4/\/3. Note that each / here denotes a two- 
dimensional integral, and that u and v depend on pi and p 2 . 
Formula (6) simplifies to 



Pr( £l2 =l) = yJJ (i uv ) n ~ 2 + (! " UV T~ 2 ) d P* d K 

T T 



(7) 



Let a, 6, and c be the vertices of T. We can partition the 
domain of integration of (7), namely the set T x T of all pairs 
(Pi > P2 ) in the triangle, into three subdomains A,B,C according 
to which vertex of T is also a vertex of Ti 2 . Because of symme- 
try, the integral (7) will be the same on each subdomain, and 
we can write 



Pr( £l2 =l) = 16 J J ((uv) n ~ 2 + (1 - uv) n ~ 2 ) dpx dp 2 



(8) 



To compute (8) we will perform a change of variables, re- 
placing pi and pi by the distances u and v. This requires that 
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we rewrite the differential dpi dp2 in terms of u, v, du and dv. In 
other words, we must compute the (4-dimensional) measure of 
the set of point pairs p, q in T x T such that the line pq cuts the 
side ab at u' between u and u + du, and cuts ac at v' between v 
and v + dv. Instead of computing this directly, we will measure 
the set A'(u,v) of pairs p,q for which u' < u and v' < v. See 
figure 39. We then obtain the desired answer by diferentiating 
the measure |A'(u, v)\ of this set with respect to u and v. 

What is the set A'(u,v)7 In order for u' < u and v' < v, 
the points p, q must be inside T\i- Furthermore for each choice 
of p, the point q must be in regions I or IV of figure 40. The 
measure |A'(Ti,t>)| is therefore the area of T12 times the average 
area of these two regions when p\ ranges over the whole T\z. If 
T\i were an equilateral triangle, then from symmetry alone it 
would follow that the average area of each region I- VI is exactly 
one sixth of the area of Tu . However, since we can make T12 
equilateral by an affine map, and every afnne map preserves the 
ratio of areas, we conclude that the same is true no matter what 
is the shape of T12. Therefore, we have 



\a\u,v)\ = |r 12 |.-|r„| = 



V? V 2 

~L6~ 



(9) 



Differentiation of the formula above with respect to u and v 
gives dpi dp2 — \uvdudv. Substituting this into (8) we get 



1 1 



Pr( £l2 = l) = 16 y J {{uv) n ~ 2 + (1 - uv) n ~ 2 ) ^ du dv 

0 0 

1 1 

= 4 J J ((uc)"- 1 + (1 - uv) n ~ 2 uv) du dv (10) 
0 0 

The evaluation of (10) presents no great problems. For the 
term (uv) n ~ l we have 

11 11 
J J{uv) n ~ x du dv = ^ J u* 1 " 1 ^ J v n - x dt^ = ( n ) 

0 0 0 0 

As for the second term (1 - uv) n ~ 2 uv, we first integrate it 
with respect to u by substituting z for 1 — uv: 

1 

r 1-2 uvdu 

0 



J(l- uv) n 

1-v 




Figure 39. 




Figure 40. 



dz 
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l 

V 



„n-l 



n - 1 



n 



l-v 



l-jl-v) 11 - 1 _ 1_- (1 - v) 
n(n — l)v 



n-l 



n 



(12) 



Now we integrate expression (12) with respect to v, substi- 
tuting w for (1 — v): 



//(!-«)• 

0 0 i 

? / I -(l-t;)"- 1 _ l-^-p)"- 1 ^ 
7 \ ra(n — l)u n / 



xn 2 uvdudv 
l 



-~y U(»-i)(i- 



n-l 



) 



l)(l-w) n 

/ / l + w + w 2 + --- + w n ~ 2 _ w n ~ x \ 

~ J \ n(n - 1) n / 

o 

i Ail ^ l \ _ j_ 

" n(n- 1) V 2 + 3 + '*' + n- l) n 2 



dw 



H n -1 1 

n(n — 1) n 2 
From (11) and (13) we get 

Pr(e 12 =l) = Yl Pr(ey=l) 

l<t<j<n 



(13) 



Vn 2 n(n-l) n 2 / 



n(n — 1) 



(14) 



We conclude that 
h(n) 
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Other distributions. Renyi and Sulanke were also able to 
compute the asymptotic expansion of h(n) for a few other dis- 
tributions. If the points are uniformly distributed over a fixed 
convex polygon K with vertices 01,03, ... ,a r , then we have 

Mn) = |(lnn + 7 )+H( £ lnj|})+o(l) (15) 

where Ti is the triangle a,_iaja, + i. 

This result may seem paradoxical: let K be a triangle, and 
suppose we remove an infinitesimally small piece K' from one 
of its corners. According to formulas (1) and (15), the average 
number of vertices in the convex hull jumps from 2 In n + o( 1 ) to 
I In n + o(l). What happens is that for small n the o(l) terms 
are still large enough to obliterate the difference between 2 In n 
and I Inn. When n is large enough to make the o(l) terms 
negligible, the probability of finding a point in K' has become 
significant — i.e., K' is no longer infinitesimally small. 

If the points are uniformly distributed on a circle, the av- 
erage value of h is 

h(n) = 2T(5/3)i^-j ^(l-Ml)) (16) 

In fact, we have h(n) — 0(n 1 / 3 ) whenever the points are uni- 
formly distributed over an arbitrary convex figure with smooth 
boundary; the shape of the curve changes only the proportion- 
ality factor. Renyi and Sulanke also show that 

h(n) = 2V27rlnn(l + o(n)) (17) 

if the points come from a two-dimensional normal distribution. 



Epilogue 

In this paper we have attempted to survey a number of 
different paradigms for the design and analysis of geometric 
algorithms. We have illustrated the use of traditional algorithm 
design methods and data structuring techniques, as well as the 
use of a number of tools that are more intrinsically geometric. 
Many algorithms employ several of these tools in combination. 
We hope that in time more and more techniques from geometry 
(as a branch of mathematics) will find use in this computational 
domain, thus linking it ever more tightly with a discipline that 
is as old as human thought itself. 
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